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'Two principal approaches used in the numerical study of unsteady inviscid compressible 
flow problems—“shock-fitting** and “shock-capturing"’—are discussed briefly. This paper 
presents a new example of the latter approach. The proposed method is monotonic and 
provides second-order accuracy in regions where the concept of order is meaningful. The 
method is simple and fast running. The monotonicity of a standard finite difference 
scheme (the well-known MacCormack scheme was chosen in the paper) is achieved by, 
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20. ASSTRACT (Continued) 

, means of a local conservative smoothing operation. It eliminated nonphysical ripples, 
which usually appear near discontinuities when nonmonotonic schemes are used. The 
method is similar to one developed by Boris and Book and gives essentially the same 
results, but is much simpler and faster. The proposed method is compared to some other 
finite difference schemes. Two one-dimensional tests are used: a linear convection 
equation and a shock tube problem. Finally, the example of nonsteady three- 
dimensional scattering of a shock wave by a step obstacle of finite width is discussed., 






CONTENTS 


0. Introduction . 1 

1. Construction of Second Order Monotonic Difference Schemes . 6 

2. Solution of the Linear Transport Equation. 11 

3. The Evolution of a Gas-Dynamic Shock . 14 

4. Calculation of Three Dimensional Interaction of a Shock Wave with an Obstacle. 19 

References . 33 












A CLASS OF MONOTONIC SHOCK-CAPTURING DIFFERENCE SCHEMES 


0. Introduction 

A broad class of gas-dynamic problems can be treated within the 
framework of ideal gas theory, i.e., assuming the gas is inviscid and 
thermally nonconductive. The differential equations of gas dynamics 
follow from the laws of conservation of mass, momentum, and energy 
in integral form assuming continuous differentiability of the fluid 
variables, and constitute a system of quasi-linear hyperbolic equations. 

It is well known that in general these may posess discontinuous solutions 
even when smooth initial data are prescribed^. Physically, the presence 
of discontinuities in the solution typically signals the appearance of 
a shock wave. There are two possibilities for correctly describing 
discontinuous flows within the framework of ideal gas theory. The 
first consists of breaking up all the regions In the problem into 
subregions of smooth flow, which are described by the differential 
equations of gas dynamics, while the discontinuities (boundaries of the 
subregions), are described by conservation conditions. It is important 
to distinguish weak discontinuities (discontinuities in derivatives of 
the fluid variables), tangential discontinuties, and finally, strong 
discontinuities (shock waves). The second approach consists in utilizing 

the conservation conditions in integral form, which allows discontinuous 
solutions. For historical reasons (the differential equations were 
derived earlier), these are referred to as "generalized," iu contrast 
to the continuous classical equations. Ve note that consideration of 

generalized solutions extends the class of possible solutions, and it 
is therefore necessary to make use of additional considerations in order 

to determine the correct solution. For example the requirement that 
entropy not decrease at a discontinuity permits us to exclude 
Manuscript submitted May 12,1981. 
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rarefaction shocks in the flow of a perfect gas (from a formal point 
of view such a discontinuity is unstable). 

In numerical modelling of gas motions in the presence of shock waves 
both of the above approaches are possible. In the former, the physical 
jump conditions on the discontinuities dividing the regions of continuous 
flow are used both in order to obtain the conditions connecting the fluid 
variables on the two sides of the discontinuity and to determine the 
motion of the computational mesh points following the discontinuity. 
Changes in the topology of the lines of discontinuity can be followed 
either using a moving curvilinear mesh or by means of an algorithm baaed 
on a fixed mesh. The first method is most completely worked out in the 
2 

paper by S. K. Godunov et al. while the second has been developed by 

3.A 

Moretti et al. 

The shortcoming of both approaches is the inhomogeneity of the 
resulting difference schemes and consequent complicated structure of lh>' 
numerical algorithm. Both approaches, as a rule, demand an a p riori know- 
ledge of the flow pattern during the calculation. 

When it is impossible to know in advance the flow pattern t or when 
it is changing qualitatively with time, it is more convenient to make t>s« 
of "shock capturing" difference schemes, amounting to difference approxi¬ 
mations to the integral form of the conservation laws for every computa¬ 
tional cell. The form of the difference equations does not depend on 
the character of the flow or the position of the possible shocks, and 
therefore such schemes are homogeneous.^’"* [Apparently "homogeneous" 
means universal, i.e., not problem-dependent.—DLB) In this approach 
shocks appear as regions of abrupt variation of the fluid quantities. 




Shock-capturing schemes can also be based on the differential equations. 
For this purpose small corrections are introduced in the equations, 
typically of nonlinear form, analogous to physical viscosity. The 
equations assume parabolic form and permit a smooth solution which 
tends to the solution of the original system as the "artificial" vis¬ 
cosity vanishes. However, difference analogs of the conservation laws 
are satisfied, in this technique, generally speaking, only approximately. 
In difference schemes based on the integral form, the mass, momentum, 
and energy conservation laws are sati^'-ed for every computational cell 
to roundoff. Such schemes are describable in terms of fluxes across 
boundaries of the computational cells, and so the conservation laws 
for each computational region are algebraic consequences of the con¬ 
servation laws for cells, i.e., the schemes are conservative.^ 

We note that the first approach, as a rule, distinguishes only 
simple shocks, while the other types of discontinuity are calculated 
as in the shock-capturing method. A difference scheme employed for 
solving practical problems is required to be accurate, that is, it 
must be able to describe flows on comparatively coarse meashes; it must 
be economical; and it should be simple to employ. The accuracy of the 
scheme for calculating smooth flows is determined by the order of the 
approximation. In the neighborhood of a shock, the change in the fluid 
variables is comparable with their magnitude, so that the concept of 
the order of the approximation becomes meaningless. It is well known 
that schemes of the first type lead to smearing of the shock because 
of the strong numerical viscosity. Schemes of second and higher order 
give rise to significantly less smearing of shock runs, but are typically 
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nonmonotonic. One of the principle reasons for the appearance of un¬ 
physical oscillations in the neighborhood of a shock is the nonvanishlng 
dispersion of the difference scheme, which leads to misrepresentations 
of the form and speed of propagation of physical disturbances, especially 
at short wavelengths. Nonmonotonicity may also be produced by the non¬ 
linearity inherent in real problems.^ 

Artificial viscosity is introduced to suppress unphysical oscil¬ 
lations. Host methods for achieving this, however, do not yield a 
monotonic scheme, and consequently make the localization of shocks worse 
and strongly smear large gradients. This problem becomes especially 
severe in calculations involving strongly shocked flow arising from 
nonsteady interaction of shock waves with one another and with obstacles. 2 
Recently a number of methods using nonlinear limiters of various sorts 

7-14 

to filter unphysical oscillations have appeared. Many of these are 

highly effective, so that the assertion that "good schemes of first and 
of higher order smear shock fronts in practically the same degree" 2 
cannot be accepted as correct. 

It is convenient to choose an actual difference scheme in several 
stages. First the dispersion and dissipation properties are studied 
in the linear approximation using the method of differential approxi¬ 
mations or harmonic analysis for the simplest transport equations. It 
is essential to test the schemes selected through linear analysis on 
simple one-dimensional problems having an exact solution. The final 
stage of development of any difference scheme must be its verification 


on model two-dimensional problems. 









In the present paper a nonlinear conservative smoothing procedure 
is proposed in order to achieve monotonicity in explicit schemes of 
higher order.' Although this smoothing procedure is applicable to 
arbitrary schemes, the underlying transport algorithm is taken here to 
be second-order MacCormack differencing.^ The results of the linear 
analysis of the dispersion and dissipation properties of this and a 
number of other widely used schemes are presented in Ref. 16. Below, 

‘ the second stage in the choice of an optimal scheme is considered in 

detail. This is conveniently divided into two parts. The first is an 
investigation of the solution of the linear advection equation, which 
^ facilitates an understanding of the nonlinear properties of the scheme. 

The second part is a calculation of the flow resulting from the evolution 
. • of a one-dimensional gas dynamic shock. In conclusion, as an illus- 

! tration of the possibilities of the proposed numerical method, a 

three-dimensional problem involving the interaction of a plane shock 
wave with an obstacle is discussed. 


1. Construction of second order monotonic difference schemes 

An effective device for constructing monotonic difference schemes 

having second-order accuracy for smooth flows was presented by Boris 
11 12 

and Book. ’ It consists of two stages. In the first a large numeri¬ 
cal diffusion is Introduced into the solution, which guarantees the 
monotonicity of the scheme. In the second stage this diffusion is 
canceled wherever doing so does not introduce new unphysical extrema or 
accentuate existing ones. By means of this approach a family 
of FCT (Flux-Corrected Transport) methods was constructed. 

We consider the simple example of the Cauchy problem for the 
linear advection equation 


where 


3f 

at 


3f 

3x 


0 


( 1 - 1 ) 


f(0,x) = f Q (x). 

The exact solution has the form f(t,x) = fg(x ” vt )• Let ^ the 
value of the function at the i**! grid point of a uniform mesh on the nth 
time level; let i + T be the operator which propagates the solution 
from the nth to the (n+l)st level. The form of T depends on the 
particular difference scheme used. We define a diffusion operator D by 


Df i ' * 1 + 1/2 - * 1 - 1/2 ' Qt£ 1 + 1/2 - Q { f i - l /2 a ‘ 2) 

where 

6f i+l/2 = f i+l " f r 
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Then the first (diffusion) stage can be represented in the form 


? 1 ■ (1+T+D)f°. 

Now we cancel the diffusion in such a way that in the resulting solu¬ 
tion no new extrema in comparison with f^ appear, and those already 
present are not amplified. For this purpose in the antidiffusion 
operator A, given by 

"i ■- ( *i+l/2 - <-!«>> 

we limit the fluxes <p* + £/2 = ^1+1/2 accordin 8 to the formula 

*1+1/2 ‘ S • “ x l0 - Bd “ lS 5? i-l/2' l*l+l/21' 

S ^ 1+3/2 l). (1.3) 

where S = sign cpj+1/2- 

The transition from the nth to the (n+l)st time level has the form 

f” +1 = (l+A)^ = (1+A)(l+T+D)f° . (1.4) 

This technique—the explicit cancellation of the diffusion—leads to 
retention of some of the diffusion in smooth regions [i.e», where the 
limiter (1.3) does not operate and A =? - D], even when T = 0. Conse¬ 
quently Boris and Book proposed two other algorithms: 
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'phoenical' 


fj +1 - I(1+A)(1+T) + D]fJ (1.5) 

and 

implicit 

f“ +1 = (1+A)~ 1 (l+T+D)f° > (1.6) 

which permit complete cancellation of the diffusion when T = 0 in 

e 

smooth regions. We note that using the algorithm (1.4)-(1.6) leads to 
substantial improvement in the results in comparison with widely used 
schemes. 

In the SHASTA* 1 scheme the coefficient Q *» 1/8. Subsequently 
Boris and Book investigated In detail the influence of the value of 
magnitude of Q on the dissipative properties of this scheme and deter¬ 
mined the optimum dependence of Q on the Courant number. 1 ^ However, 
practically speaking, this optimization is meaningful only for a linear 
equation with constant coefficients. For the equations of gas-dynamics 
as calculations reveal, good results are obtained for Q between 
1/10 and 1/6. 

One can point out three factors responsible for the success of the 
method of Boris and Book. First, the diffusion stage, as shown by 
linear analysis, substantially improves the dispersion properties of 
the scheme. Second, the diffusion is introduced in a conservative 
manner, and, finally, it is done so as to permit excellent localization 
in the neighborhood of the shocks. Evidently it is precisely these 
last two properties which are most important. In Ref. 13 it was . 



proposed to interpret the method of Boris and Book as a nonlinear con¬ 
servative smoothing procedure and to introduce diffusion at the (n+l)st 
time level according to 


(1+A)(1+0)(1+T)fJ 

(1.7) 

(1+A+D)(l+T)f° 

(1.8) 


This procedure fails to improve the dispersive properties and falls 
somewhat short in quality of the results of the algorithms (1.4)-(1.6). 
However, in contrast to the original method of Boris and Book, which is 
essentially one-dimensional in character and which in two- 
dimensional calculations can only be applied via the method of coordi¬ 
nate timestep-splitting of the the difference equations, smoothing in 
the form of Eqs. (1.7)-(1.8) is easily incorporated in essentially 
two-dimensional schemes. For this purpose, after each timestep successive 
smoothing operations are carried out with respect to each of the coordi¬ 
nates. Such smoothing has turned out to be effective in solving a 
whole series of problems.* 7 

In the present paper a simpler smoothing procedure is proposed, 
having a local character while retaining the conservative property of 
the difference scheme. A constant diffusion is introduced only in 
regions which are not monotonic, i.e., where the solution has a numerical 
ripple. As a result of this, just as in the method of Boris and Book, 
there is a certain amount of smearing of physical extrema, but for the 
most part only oscillations produced by the numerical scheme are 










removed. The conservative property is assured by introducing the 
diffusion in the form of fluxes across the cell boundaries. Using 
the notation introduced above, ve can write the smoothing at the nth 
time level as 


f" +1 - (1+T)fj + D*fJ, (1.9) 

and on the (n+l)st time level in the form 

l 

f^ 1 - (l-HJ*)(l+T)f^. (1.10) 

In these expressions 

"X ' 4/! - il/2' 

where 

* ( ^i+1/2’ 6f i+l/2 6f i+3/2 < 0 or 6f i+l/2 6f i-l/2 < 0 

Vl/2- (1 * n ? 

I 0 otherwise. 

In order to code this prescription it is convenient to search the entire 
array frecording the boundaries across which the flux is nonvanishing, 
and then in a second pass calculate the nonzero fluxes. 

This procedure for introducing artificial viscosity runs substan¬ 
tially faster than the method of Boris and Book, while yielding com¬ 
parable results. Note that, in contrast with the algorithms (1.4)-(1.7), 
smoothing according to the prescription (1?8)-(1.10) in regions where 
the fluid quantities vary monotonically in general leaves the solution 
unchanged. 
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2. Solution of the linear transport equation 

In this and the following section the proposed numerical method 

is compared with a number of different schemes, using a one dimensional 

22 

model problem having an exact solution. The schemes of S. T. Godunov , 

15 23 

MacCormack , V. V. Eremin and Yu. M. Lipnitskii , of first, second 

and third order, respectively, are considered and also the non-linear 

scheme of Boris and Book, SHASTX 12 , embodying the algorithm (1.5). These 

schemes are widely employed and have performed well in practice; therefore 

the present analysis permits one to carry out an objective evaluation of 

the proposed method. 

Let us approximate equation (1.1) on a uniform mesh with Lx. m 1. 

The computational mesh contains a total of 100 cells. On the boundaries 
periodicity conditions are imposed. The initial conditions are given in 
the form 

( 0.5, x > L 

<p00, lsxsL 

where L = 20 and 

cp(x) ■ 2 (2.1) 

or 

<p(x) = 0.5 + 0.075x. (2.2) 

In Fig.. 1 are shown the results of the calculations after 800 
time steps with Courant number v= 0.2, i.e. for a total running 
time t - Vt/L, after which the initial square wave (2.1) of width 20 . 
cells has been carried across 160 cells. The curves (al), (a2), (a3), 
correspond to the schemes of Godunov (which coincides in the linear 
case with one-sided differencing), MacCormack (equivalent to the scheme 
of Lax and Wendroff), and Eremin-Lipnitskii; (b). to the SHASTX scheme; (c). 
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(d), (e) correspond Co Che monotonic MacCormack scheme with smoothing 
according to algorithms (1.7), (1.10), and (1.9). The best results 
were obtained with SHASTX. The SHASTA* -1 algorithm and the monotonic 
MacCormack scheme with the smoothing (1.9) were very slightly inferior 
to this. The smoothing procedures (1.7) and (1.10) give worse results, 
which are comparable, however, with the results obtained from the third- 
order scheme. The first-order and second-order schemes yielded 
unsatisfactory approximations to the solution. 

Since the solutions of the linear advection,equation using the 
algorithms (1.7) and (1.10) were inferior, however slightly, to the algorithm 
(1.9), in the remainder of this section we will present the results of 
calculations carried out using the latter. We must point out, however, 
that these other algorithms will be analyzed below, inasmuch as they 
possess definite advantages in solving real gas-dynamic problems. The 
results of the present section are applicable in problems where linear 
equations are solved, for example, in certain meteorological problems. 

Figure 2 illustrates the "flexibility” of the various schemes in the 
sense of Ref. 2, i.e., the capability of solving problems in which the 
Courant number varies greatly over the computational region. Here 
profiles are shown at t * 3, obtained for Courant number v» 0.6 
(curve 1) and v= 0.2 (curve 2) using first-order (a), second-order (b), 
third-order (c) and the monotonic MacCormack (d) schemes. The 
superiority of the latter is indubitable. The calculation was not 
carried out with SHASTX, which requires that condition v< 0.5 be 
satisfied 11 . 

As we remarked in Section 1, the smoothing which removes numerical 





oscillations smears out physical extrema. Consequently, the test 
problem using initial conditions in the form of a triangular profile 
(2.2) constitutes a stringent test of a scheme. Figure 3 shows 
profiles at time t * 1 (curves 1) and t ■ 20 (curves 2) for the same 
schemes as in Fig. 2, if v- 0.2. The first-order scheme practically 
"forgets" the Initial conditions; the second-order scheme qualitatively 
changes the shape of the pulse as time elapses, while the thid-order 
scheme preserves the "height" of the dlstrubance well. The monotonic 
MacCormack scheme is only negligibly inferior to SHASTX (the dashed lines 
in Fig. 3d). 
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3. The evolution of a gas-dynamic shock 

The problem of the evolution of a gas dynamic shock Involves all the 
characteristic features of one-dimenrional ideal gas flow: creation and 
propagation of a shock wave, contact discontinuity and rarefaction fan. 

By comparing the results of the calculation with the analytic solution, 
we can check both the accuracy with which the characteristic flow regions 
are transported and the speed with which they develop-. A comparison 

of the various schemes mentioned above will be presented, based on this 

) 

problem. 

As was stated earlier, unphysicaloscillations in the neighborhood of 
shocks and steep gradients appear because of the nonvanishing dispersion 

associated with a difference scheme, the nonlinearity of the gas-dynamic 

24 

equations, and also for several other reasons . Besides this, it is 
possible that the various nonlinear smoothing operators may Interact 
with one another nonselfconsistently. 


The system of equations has the following form: 




+ P) 



C 2 \ 2 C 2 f 2 

J e I dx + J (e + p )u I < 


Here e - p/(y-l) + pu /2, where y is the adiabatic index. We consider 


the evolution of a shock with initial pressure ratio p^/p^ “ 2 and 
density ratio pj/p^ " an< * Y “ 1.4. In Fig. 4 are shown density 
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profiles obtained using the Godunov (a), Eremin-Lipnitskii (b), 

MacCormack (c) and SHASTX (d) schemes. The exact solution is 
indicated by the heavy continuous line and the numerical results 
up to t s 25 by a fine continuous line. The first-order scheme, as 
one might expect, badly smears the shock front and contact discontinuity 
and represents the rarefaction fan poorly. The Eremin-Lipnitslcii and 
MacCormack schemes are nonmonotonic and yield oscillations in the 
neighborhood of the shocks. SHASTX, which represents the shock front 
very well, gives rise to small undershoots in the density in the neighborhood 
of the contact discontinuity and rarefaction fan. For this case the 
SHASTA scheme, as described by Boris and Book, gives practically the 
same results as does SHASTX. All of the schemes accurately represent 
the pressure profile between the shock wave and the rarefaction fan. 

In Fig. 5 are shown results of calculations using the monotonic 
MacCormack scheme with various forms of smoothing. Fig. 5a corresponds 
to algorithm (1.7). One shortcoming of this approach lies in the presence 
of undershoots in density in the neighborhood of the shocks and rarefaction 
fans. Very similar results are obtained when one applies the smoothing 
(1.9) to all of the equations in system (3.1) (Fig. 5b). 

In the present work it is proposed to carry out self-consistent 
smoothing of the density, momentum, and energy in regions where 
nonmonotonicity in the density profile exists, using the procedures 
(1.9) and (1.10) described above. The justification for this lies in 
the fact that every stable one-dimensional shock has a density 
discontinuity. In particular, this self-consistent smoothing 
introduces no changes in the velocity at a contact discontinuity, since 







the density and radmantum are smoothed in a "related" manner. Zt also 
somwehat reduces the expenditure of machine time, since in order to 
determine nonmonotonicity, one array of numbers instead of three is 
searched [Cf.(l.ll)j. We note that this type of smoothing cannot be 
implimented using algorithms (1.4) - (1.8) because of their two-stage 
character. 

The density profiles obtained using self-consistent smoothing at 
the n th and (rri-1) st levels are shown in Figs. 5c and 5d. These 
results are clearly superior to all those shown previously. It is 
interesting to observe that the difference in the results obtained by 
smoothing according to the algorithms (1.9) and 1.10) is much smaller 
than in the case of linear advection (Fig. 1). This is indicative of 
the decisive role played by nonlinear effects in numerical solutions 
of the gas-dynamic equations. 

We pause to draw attention to some peculiarities in the solutions 
we have found. Most of the calculations were run with At " 0.5. 

SHASTA, SHASTX, and the MacCormack scheme with self-consistent 
smoothing, however, yielded a highly oscillatory solution. The 
oscillations are a consequence of the nonselfconsistency of the smoothing 
and the interaction between the resulting errors in the various fluid 
quantities. A separate paper will be. devoted to a detailed analysis 
of the nonlinear properties of the different schemes, the nonmonotonic*ty 
resulting from various causes, and ways of contending with the unphysical 
oscillations that result. We therefore will not go inot the reasons for 
the appearance of computational ripples in detail, but content ourselves 
with remarking that they vanish as the time step is reduced. The • 
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solutions shown in Fig. 4c and 5b were obtained with a time step 
t - 0.25. 

We now compare the respective running times. Taking as the unit of 
time that required for a single step according to the HacCormack scheme, 
the corresponding running times for the other schemes ate as follows: 
monotonic MacCormack with self consistent smoothing, 1.4; with 
nonselfconsistent smoothing according to algorithm (1.9), 1.5; using 
algorithm (1.7), 2; for the Godunov scheme 2 - 2.5, depending on how 
the "large" quantities are calculated; for SHASTA and SH/.'.TX, 3; and 
for the Eremln-Lipnitskii scheme, 2.5. 

We next consider an example in which a strong shock wave is formed. 
The initial parameters are taken to he the same as in the paper of Boris 
and Book 11 (p 2 /p 1 = 480, ” 8 » Y “ 5/3, At. « 0.02). The Erernin- 

Lipnitskii and MacCormack schemes in the absence of an additional 
artificial viscosity do not permit calculation of such a flow. SBAvr/. 
and the monotonic MacCormack scheme with smoothing according to 
algorithm (1.7) give approximately the same results as SHaSTX. We 
therefore make further comparisons of the results obtained, using the 
method of Godunov, SHASTX, and menotonic MacCormack scheme with h> f- 
consistent smoothing given by (1.9). In Figs, 6 and 7 a) c shown density 
and pressure profiles obtained using these schemes at t.iw- 1^2 (deshf-d 
lines and t = 4 (fine continuous lines). A vertical bashed-dotted line 
is used to show the initial position of the shock. The virtues and 
shortcomings of the numerical methods under invest igat ion rare 
clearly shown i*. Table 1, in which the values of density and relative 
errors at time t ■ 4 are shown. In the left column of the table is shown 
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the number of computational cells. The Godunov scheme strongly 
smears the shock front and contact discontinuity and approximates the 
exact solution poorly in the region of the rarefaction fan. SHASTX 
is somewhat better than the monotonic MacCormack scheme in representing 
the shock wave and rarefaction fan, but somewhat inferior to it in 
representing the contact discontinuity. In addition, SHASTX gives 
incorrect values in the neithborhood of the initial shock. 

The monotonic MacCormack scheme and SHASTX have good dynamical 
properties. They permit relatively fast determination of the true flow 
pattern (Cf. the Godunov scheme). This last property is especially important 
in solving nonstationary problems. 
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4. Calculation of three dimensional interaction of a shock wave with 
an obstacle. 

The possibilities of the proposed method are illustrated below for 
the example of reflection of a plane shock wave from a obstacle of 
finite width. The initial stage of the interaction precess is 
investigated, in which the reflected shock wave begins to form and the 
pressure loading on the object is a maximum. We remark that this problem 
is of interest in its own right, inasmuch as experimental determination of 
nonstationary three-dimensional flows encounters serious difficulties. 

The problem is solved in Cartesian coordinates. The mesh is 
uniform. Timestep splitting was employed in solving the difference 
equations 2 ^* 2 ** which permitted a substantial saving in running time. 

The conservation laws in Eulerian coordinates have the form 


i ds dt e 0, 


j P | 2 dv + J t2 pV-n 
v C 1 fc l s 

•r* (pw + p) 

v t 1 fc l S 

/• l; dv + a (e + p) V • n ds dt 


n ds dt = 0, 


We denote by L q (At) the operator which performs a step in the direction o: 
n + 1/3 


ijk 


• L (At) £ 


ijk. 


L is defined in the following manner: 
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where 


T ljk • f ljk - f ^s- r ijk). 


n+1/3 „ I r f n + 7 _ M. f p )] 

r ijk 2 L ijk r ijk Ah u ijk i-q,j-r,k-s' J » 


q"5,r«*6,s = 6 

xa ya za. 

The smoothing is carried out according to the algorithm (1.9). The 
transverse momentum components were smoothed independently. 

The complete transport operator will look as follows: 


Written thus, even if the split operators L o (At) not commute, the 

27 

difference scheme has second-order accuracy over a double timestep 
The stability condition for the split equations is 
At £min(At ), o * x,y,z. 


At • £ min (—t~-) » 

° ijk V* 

and a is the speed of sound. This condition turns put to he 
considerably less stringent than for the unsplit operators: 


At & min ( ——- ;-) » 

Ivl + a/3 
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Since At^ ( x is the direction of propagation of the shock wave) 
is usually much smaller than At^ and At z , it makes sense to carry out 
the calculation in the y and z directions with a double tlmestep. Thus 
using 

L a (2At)~ L a (At) L (At), a = y,z, 
we finally obtain 

n+2mAt . 

f = L (At)[L (At) L (2At) L (At) L (2At)]® -1 

ijk y x z x y 

. L (At)L (2At)L (At)L (At)f" . 
x z x y lj k 

Use of this algorithm instead of 

^ n+2mAt n 

f = L®(2At)f 
ijk 

permits reduction of the computational volume by a factor of 1.5. 

In the test the computationalmesh consisted of 20 X 20 X 15 cells. 

Calculation of a single case up to t = 1.2 - 1.5 required about 2 hours of 

time on the M4030. Here t = tV /h, where V is the speed of the 
amp amp 

reflected shock wave which arose as a result of the interaction of the 
incident wave with the square object,, and h is the step height. 

A typical flow pattern is shown in Fig. 8. At t = 1.1 the isobars 
are shown in (a) with separation Ap = 10.0,and the lines of constant 
M*ch number (b) at intervals AM = 0.2. The Mach number of the incident 
shock wave =5.0 and the ratio of width to height of the step was 
f/h = 5.2. The adiabatic index was y = 1.4 and tha pressure in the 
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undisturbed gas was = 1.0. It is clear that the method permits 
excellent resolution of the shock wave and the flow in the neighborhood 
of the corner points. 

In Fig. 9 are shown the isobars in two cross sections x * const, 
located equidistantly and parallel to the boundary of the step at a 
distance 0.3 h. The ratio I/h - 2.0 in (a), 3.6 in (b), and 5.2 in (c). 

In Fig. 9a are shown the Isobars for the plane case (If h - »). 

Figure 10 illustrates the influence of the ratio i/h on the bending 
of the shock wave in the xy plane. Curves 1,2,3 ^correspond respectively 
to 1/h = 2/0, 3.6, and 5.2 (the bending is related to the bending for the 
plane case). 

The authors thank Yu. P. Golovachev for his unstinting attention to 
this work and useful discussions, and also V.M. Goloviznin, G.L. Stenchikov, 
and A.P. Favorskii for valuable suggestions. 
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Olcy Attn R Plebuch 

TF.W Defense 6 Space Svs Group 
P 0 Bo: 1310 

San Bc-rnardino, CA 92402 
Olcv Attn E Wong 
Olcy Attn P Dai 

Universal Analytics, Inc. 

7740 W Manchester Blvd 
Plava Del Rev, CA 90291 
Olcy Attn E Field 

Weidlinger Assoc., Consulting Eng 
110 E 59th Street 
New York, NY 10022 
Olcy Attn M Baron 

Weidlinger Assoc., Consulting Eng 
3000 Sand Hill Road 
Menlo Park, CA 94025 

Olcy Attn J Isenberg 

Westinghouse Electric Corp. 

Marine Division 
Kendy Avenue 
Sunnyvale, CA 94088 
Clcy Attn W Volz 
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